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Abstract 

In this paper we study the depth at which a cosmic ray shower 
reaches its maximum (X mSLX ) as predicted by Monte Carlo simulation. 
The use of J max in the determination of the primary particle mass 
can only be done by comparing the measured values with simulation 
predictions. For this reason it is important to study the differences 
between the available simulation models. We have done a study of the 
first and second moments of the A max distribution using the Corsika 
and Conex programs. The study was done with high statistics in the 
energy range from 10 17 to 10 20 ' 4 eV. We focus our analysis in the differ- 
ent implementations of the hadronic interaction models Sibyll2.1 and 
QGSJetII in Corsika and Conex. We show that the predictions 
of the (X max ) and RMS(X 

max) depend slightly on the combination of 
simulation program and hadronic interaction model. Although these 
differences are small, they are not negligible in some cases (up to 5 
g/cm 2 for the worse case) and they should be considered as a system- 
atic uncertainty of the model predictions for (A max ) and RMS(A max ). 
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We have included a table with the suggested systematic uncertainties 
for the model predictions. Finally, we present a parametrization of 
the X max distribution as a function of mass and energy according to 
the models Sibyll2.1 and QGSJetII, and showed an example of its 
application to obtain the predicted X mSLX distributions from cosmic 
ray propagation models. 

1 Introduction 

The cosmic ray composition at the highest energies is probably the most 
difficult and most meaningful question yet to be solved in the present as- 
troparticle physics scenario. Due to the unknown strength and structures of 
the magnetic fields in the Universe, anisotropy studies are also intrinsically 
dependent on the mass composition and a better identification of the sources 
is probably only going to be possible if the cosmic ray composition is known 
beforehand. 

The most reliable technique to infer the mass composition of showers with 
energy above 10 17 eV is the determination of the X max and posterior compar- 
ison of the measured values with predictions from Monte Carlo simulation. 
This is because above 10 17 eV fluorescence detectors can measure X max with 
a resolution of 20 g/cm 2 . The evolution of the detectors, the techniques 
used to measure the atmosphere, the advances in the understanding of the 
fluorescence emission and the development of innovative analysis procedures 
have resulted in a high precision measurement of X max and its moments. 
The Pierre Auger Observatory [?], the HiRes Experiment [?], the Telescope 
Array [?] and the Yakutsk array [?] quote systematic uncertainties in the 
determination of the (X max ) to be 12, 3.3, 15 and 20 g/cm 2 , respectively. 
Considering the quoted errors and taking into account that the data have to 
be compared to simulation predictions, it is very important to understand 
the details and reduce the differences between simulation programs. The 
proposed experiment JEM-EUSO [?] is also going to use the fluorescence 
technique to detect air shower from the space. For this reason, we extended 
the analysis done in this work up to 10 20 ' 4 eV within the energy range aimed 
by JEM-EUSO. 

The dependency of X max with primary energy and mass (A) has been 
analytically studied in a hadronic cascade model [?]. Monte Carlo programs 
can simulate the hadronic cascade in the atmosphere using extrapolation 
from the measured hadronic cross sections at somewhat lower energy. It has 
been shown before that different hadronic interaction models do not agree in 
the prediction of the (X max ) and other parameters [?]. 
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In this paper we study in detail the dependence of (X max ) and RMS(X max ) 
as a function of energy and primary mass. We compare the result of two 
hadronic interaction models. We have done a high statistics study and we 
show that the discrepancies between models and programs are at the same 
level of quoted systematic uncertainties of the experiments. The analysis 
done here points to the need of a better understanding of the interaction 
properties at the highest energies which can be achieved by ongoing analysis 
of the LHC data which already resulted in updates of the hadronic interac- 
tion models. At the same time the results presented in this paper point to 
discrepancies between different implementations of the same hadronic inter- 
action model which need to be better understood. 

We also present a parametrization of the X max distribution as a func- 
tion of mass and energy. Several theoretical models have predicted the mass 
abundance based on astrophysical arguments [?, ?, ?, ?]. In order to compare 
the predicted abundance with measurements, one has to convert the calcu- 
lated flux for each particle into X mauX . Until now, this could only be done 
using full Monte Carlos simulations. We present here a parametrization of 
the X max distribution to allow the conversion of astrophysical models into 
A max measurements. Parametrizations of (X max ) as a function of energy and 
mass have been already studied [?]. What we present here is a step forward, 
we show the parametrization of the X max distribution which is good enough 
to calculate the first and second moments of the distribution. 

In section |2j we study the dependence of the results regarding simulation 
limitations like thinning and nucleon types. In section |3j we compare the 
models. In Section [4] we show the parametrization of the X max distribution. 
Section [5] concludes our analysis. 

2 Shower Simulation 

In this work we have used Conex [?, ?] and Corsika [?] shower simulators. 
Conex uses a one dimensional hybrid approach combining Monte Carlo sim- 
ulation and numerical solutions of cascade equations. Corsika describes the 
interactions using a full three dimensional Monte Carlo algorithm. By using 
analytical solutions, Conex saves computational time. On the other hand, 
Corsika makes use of the thinning algorithm [?, ?] to reduce simulation 
time and output size. 

Both approaches have negative and positive features. Corsika offers a 
full description of the physics mechanisms and a three dimensional propaga- 
tion of the particles in the atmosphere. However, it is very time consuming, 
limiting studies which depend on large number of events at the highest en- 
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ergies. The thinning algorithm introduces spurious fluctuations that have to 
be taken into account in the final analysis. Conex is fast, but on the other 
hand it offers only a one dimensional description of the shower. The use of 
intermediate analytical solutions might also reduce the intrinsic fluctuation 
of the shower. In the following sections both programs are compared in detail 
concerning the X max calculations. 

The hadronic interaction models for the highest energies were developed 
independently of the programs that describe the showers. For each shower 
simulator many hadronic interaction models are available. We have used 
QGSJetII.v03 [?, ?] and Sibyll2.1 [?] in this work. For the low energy 
hadronic interaction we have used GHEISHA [?] in all simulations. 

Showers have been simulated with primary energy ranging from 10 17 ' to 
10 20 ' 4 eV in steps of log 10 (E/eV) = 0.1. We have simulated seven primary 
nuclei types with mass: 1, 5, 15, 25, 35, 45 and 55. For each primary particle, 
primary energy, and hadronic interaction model combination, a set of 1000 
showers has been simulated. The zenith angle of the shower was set to 60° 
and the observation height was at sea level corresponding to a maximum 
slant depth of 2000 g/ cm 2 allowing the simulation of the entire longitudinal 
profile of the showers. The longitudinal shower profile was sampled in steps 
of 5 g/cm 2 . The energy thresholds in Corsika and Conex were set to 1, 1, 
0.001 and 0.001 GeV for hadron, muons, electrons and photons respectively. 



2.1 Fitting the longitudinal development of the shower 

For all studies in this paper, X max was calculated by fitting a Gaisser-Hillas [?] 
function to the energy deposited by the particle through the atmosphere. We 
chose a four parameter Gaisser-Hillas (GH4) function given by: 

d£. v . dE™' ( X-X y^F^ fX mm -X\ 

in which ^ max ^ \ and X max are the four fitted parameters and X is 
the slant atmospheric depth. The first guess of the X max parameter in the 
fitting procedure was chosen to be the maximum of a three degree polynomial 
interpolated within the three points in the longitudinal profiles with largest 
The full simulated profile was fitted. 

We have studied the effect of fitting a different function to the longitu- 
dinal profile. Instead of a Gaisser-Hillas function with four parameters, we 
have also fitted a Gaisser-Hillas function with 6 parameters (GH6). In this 
function, A is defined asA = axX 2 + 6xX + cin which a, b and c are also 
fitted. The differences in X max and RMS(X max ) calculations for both fitted 
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functions were smaller than 2.5 g/cm 2 and 0.6 g/cm 2 , respectively. Anyway, 
there might be a systematic effect in the determination of the X max due to the 
fitting procedure and the fitted function chosen to describe the longitudinal 
profile, which is not investigated in this paper. 

2.2 Thinning analysis 

In order to save time and output size, Corsika uses a thinning algorithm [?, 
?]. The thinning factor fthin defines the fraction of the primary particle 
energy Eq below which not all particles in the shower are followed. Particles 
with energy below E cut , where E cut = fthin * E , are sampled, some are 
discarded and others followed. Each active particle in Corsika has a weight 
attribute which compensates for the energy of the rejected ones such as that 
energy is conserved. 

The thinning algorithm causes artificial fluctuations in the calculation of 
the shower development which needs to be taken into account. For example, 
figure [T] shows the longitudinal development of one shower simulated with 
three thinning factors 10 -5 , 10~ 6 and 10 -7 . The first interaction altitude 
was fixed at 60 km and the target is Nitrogen nuclei. Sibyll2.1 was used 
for this study. It illustrates how the fluctuations of the simulated longitudinal 
profile increase with increasing thinning factor. 

Figure [2] shows the X max distribution for 1000 simulated shower for three 
thinning factors. In this example the first interaction point and target were 
not fixed. Sibyll2.1 was used again for this study. Figure [2] shows that the 
(X max ) and RMS (Xmax) are very similar for any thinning factor used. The 
maximum difference is 0.4 g/cm 2 for (X max ) and 2.8 g/cm 2 for RMS(X max ). 
Based on this study we chose to simulate all showers with thinning factor 
10" 5 . 

2.3 Isobaric Nuclei Analysis 

Some of the simulation used in this paper have been produced for another 
study and have been re-used in this work to save computational time. The 
previously simulated showers have exotic primary particle with unstable num- 
ber of protons and neutrons. In this section we compare the longitudinal 
development of showers started with exotic nuclei with the shower started 
with stable nuclei with the same total number of nucleons. Our intention is 
to show that the development of a shower at high energies does not depend 
on the number of protons and neutrons independently but depends only on 
the total number of nucleons considered. 
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Corsika and Conex simulators differentiate isobaric nuclei by allowing 
the determination of the number of protons and neutrons of the primary nu- 
clei. However the treatment of the first interaction is done by the hadronic 
interaction models. QGSJetII does not differentiate isobaric nuclei inter- 
actions, Sibyll2.1 does differentiate. 

At the highest energies the energy loss of protons and neutrons is negligi- 
ble when compared to the total energy and therefore only the total number 
of nucleons should influence the development of the shower. On the other 
hand, Coulomb dissociation should also be taken into account [?]. Neverthe- 
less, none of the hadronic interaction models available include this effect and 
therefore the development of the simulated shower should not depend on the 
number of protons and neutrons in the nuclei. 

Figure [3] shows the X max distribution for 10 19 eV showers. We have sim- 
ulated nuclei with different numbers of proton and neutron constituents. We 
can conclude that, at the energy range of interest, the predicted (X max ) and 
the RMS(X max ) does not depend on the number of protons and neutrons 
which form a nucleus with given mass A. The maximum difference in the 
(X max ) was 2.9 g/cm 2 for A = 24 and QGSJetII and in the RMS(X max ) 
was 2.5 g/cm 2 for A = 1 and Sibyll2.1. 



3 Results 



Figures 111 and [5] show the comparisons of (X max ) using both simulators for 



Sibyll2.1 and QGSJetII respectively. Figures 4e and 4f show the differ- 
ence between Corsika and Conex predictions as a function of energy and 
mass respectively when both programs used Sibyll2.1. The same is shown 
in figures 5e and 5f for QGSJetII. The differences between the (X max ) pre- 
dicted by CORSIKA and CONEX are smaller than 7 g/cm 2 in the parameter 
space studied by us. Conex tends to simulate showers slightly deeper than 
Corsika. 

Similar results are presented in figures [6] and [7] for the RMS(X max ). In the 
parameter space studied by us the differences in the RMS(X max ) calculated 
by Corsika and Conex are smaller than 8 g/cm 2 . No significant trend of 
the difference with energy or mass was seen. 

The evolution of RMS(X max ) shown in figure 6a, 6b, 7a and 7b shows 



large fluctuations apparently larger than the estimated statistical fluctua- 
tion. The statistical fluctuation shown as error bars of RMS(X max ) is the 
standard statistical variance of the variance of a distribution. No Gaussian 
approximation was used. The trend of RMS(X max ) with energy is statisti- 
cal incompatible with a linear behavior. A linear fit of RMS(X max ) versus 
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energy shows a mean x 2 /NDOF = 123/35 ~ 3.5. 

Figure |8] summarizes the differences in (X max ) and RMS(X max ) between 
the simulation programs and between the hadronic interaction models. This 
figure shows simultaneously the (X max ) and RMS(X max ) values, where the 
corresponding (X max ) and RMS(X max ) for a nuclei with mass 55 has been 
taken as reference (as suggested in [?]). This figure illustrates the importance 
of taking into account the simulation program differences into the systematic 
uncertainty of the model predictions. Each blob corresponds to the (X max ) 
and RMS(X max ) predictions for one primary particle at different energies. 

The elongation rate theorem [?, ?, ?] proposes the use of the slope of 
the variation of (X max ) with energy as a composition parameter. According 
to this proposal, changes in this slope represents changes in composition. 
Figure [9] show the slope of the variation of (X max ) and RMS(X max ) with 
energy as a function of primary particle mass for Corsika and Conex using 
Sibyll2.1 and QGSJetIL The slope in this figure corresponds to the fit of 
a straight line in the energy range (10 17 < E < 10 20 ' 4 eV). There is a good 
agreement between Conex and Corsika when the same hadronic model 
is used. The slope of the RMS(X max ) when Sibyll2.1 is used presents the 
largest discrepancy between Conex and Corsika (see|9}c). 

It has been shown before [?] that the dependencies of the (X max ) and 
RMS(X max ) are not strictly straight lines however the departure of the linear 
dependency is very small for energies above 10 17 eV as studied here. 

It is clear from figure [9] that the slope of the (X max ) as a function of 
energy is dominated by the hadronic interaction model rather than by the 
shower simulation. 



4 Parametrization of the X max distributions 



The X max distributions can be described by a function which is a convolution 
of a Gaussian with an exponential [?]: 

dX max /io-t^Mp, f t -t + a 2 /X \ 



dN 1 ^ \ X 2A 2 / V V2a 

This equation has four parameters. Nf is a normalization factor which 
gives the total number of events in the X max distribution. A, t a and a are 
parameters which are related to the decay factor of the exponential, the max- 
imum of the distribution and the width of the distribution respectively. Erfc 
is the error function. We used this equation to fit the X max distributions of 



all mass and energies we have simulated. Figure 14 shows examples X, 
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distributions we fit with this equation. The aim of this study is to use equa- 
tion^ to fit the X max distribution and calculate the (X max ) and RMS(X max ). 
We show below that a convolution of a Gaussian with an exponential allows 
a good description of the (X max ) and RMS(X max ). The proposed function is 



also a fairly good description of the X max distribution, see figures 14 



After that, we parametrized A, t and a as a function of primary mass 



and energy using the simulated showers. Figure 10 show the evolution of the 
three parameters with energy and mass. 



Figure 10 1 shows how t a has a very smooth dependence with mass and 
energy, recovering the already explored dependence of X max with mass and 
energy. On the other hand, o and A are not completely independent pa- 
rameters. Both parameters influence the width of the X max distribution in 
different ways. The parameter A changes the width of the X max distribu- 
tion by modifying the decays of the exponential, making the high X max tail 
longer or shorter, a also changes the width of the X max distribution by 
modifying the width of the central part. In fact, note that mathematically 
RMS(X max ) = v^T+A). 

Given the degeneracy in shaping the width of the X max distribution, the 
parameters a and A are inversely correlated. The parameters a and A com- 
pensate each other, fluctuations to higher values of a are correlated to fluc- 
tuations to smaller values of A. 



We performed a fit to plots in figure [10] with a linear dependence on 
logio(A) and log w (E) following equation: 

to ) 

a \ = C x X log w {E/eV) + C 2 x log 10 (A) + C 3 (3) 
^ J 

Tables [I] and [2] show the fitted parameters for Conex and Corsika 
respectively. Despite the fluctuations of a and A a linear fit in logio(A) and 
log w (E) is reasonably good approximation to describe the X max distribution. 



This can be seen in figures 11 and 12 where we show a comparison between 



the simulation, the direct fit of the X max distribution using equation [2] and 
the calculation using equation [3] and table [TJ 

It is clear that a direct fit of the X max distributions with equation [2] (blue 
lines) leads to a very good description of the first and second moments of 
the X max distribution. For all simulations and hadronic models in the entire 
energy range and for all primary particle used in our study the direct fit 
resulted in a difference on the simulation smaller than 2 g/cm 2 for (X max ) 
and smaller than 4 g/cm 2 in the RMS(X max ). 

The fit to a plane log w (A) and log w (E) is a reasonably good approxi- 



mation to describe the X max distributions. This can be seen in figures 11 
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and 12 where we show a comparison between the simulation, the direct fit of 
the X max distribution using equation [2] and the calculation using equation [3] 
and table [Q 

The parametrization as a function of energy and mass of the parameters 
that describe the Xmax distributions (equation |3| introduced some system- 
atic errors in (X max ) for the case of QGSJetII model (up to 10 g/cm 2 ). 



This is shown with red lines in figure 11 (right hand side plots). The reason 
for this systematic errors is because the parametrization used (equation [3]) is 
not the optimum one for QGSJetII model. 





Had. Model 


C x (± err) 


C 2 (± err) 


C 3 (± err) 


to 


QGSJetII 


53.06 (0.05) 


-28.74 (0.12) 


-275.93 (1.18) 


Sibyll2.1 


60.48 (0.07) 


-38.48 (0.13) 


-402.80 (1.22) 


a 


QGSJetII 


-0.26 (0.06) 


-5.63 (0.21) 


31.68 (3.38) 


Sibyll2.1 


-1.09 (0.07) 


-5.28 (0.19) 


44.41 (1.54) 


A 


QGSJetII 


-2.68 (0.14) 


-19.50 (0.43) 


100.32 (2.63) 


Sibyll2.1 


-2.61 (0.11) 


-17.89 (0.14) 


96.28 (1.76) 



Table 1: Fitted coefficients (equation |3j) — CONEX. All values in g/cm 2 





Had. Model 


Ci (± err) 


C 2 (± err) 


C 3 (± err) 


to 


QGSJetII 


53.32 (0.30) 


-29.47 (0.52) 


-283.93 (5.62) 


Sibyll2.1 


60.77 (0.23) 


-38.88 (0.31) 


-408.88 (4.67) 


a 


QGSJetII 


0.06 (0.002) 


-5.06 (0.17) 


35.99 (3.21) 


Sibyll2.1 


-0.56 (0.08) 


-4.70 (0.21) 


44.01 (2.03) 


A 


QGSJetII 


-1.73 (0.15) 


-20.63 (0.34) 


82.69 (3.54) 


Sibyll2.1 


-2.49 (0.22) 


-19.54 (0.34) 


96.04 (3.46) 



Table 2: Fitted coefficients (equation pi) - Corsika. All values in g/cm 2 . 



5 Conclusion 

We have studied the simulation programs Corsika and Conex with the 
hadronic interaction models Sibyll2.1 and QGSJetII. We have shown that 
the (X max ) and the RMS(A max ) depend slightly on the combination of pro- 
gram and hadronic interaction model chosen. It is widely known that (X max ) 
and RMS(X max ) predicted by Sibyll2.1 and QGSJetII are different mainly 
due to the different extrapolations of the hadronic interaction properties to 



9 



the highest energies. We have quantified here the differences between Cor- 
SIKA and CONEX by predicting the (X max ) and the RMS(X max ) using the 
same hadronic interaction model. These differences are small, but should be 
considered as systematic uncertainties of the model predictions. 

t ) and RMS(X n 



Figure 13 shows the evolution of the (X max ) and RMS(X max ) with en- 
ergy. No clear dependency of the difference between Corsika and Conex 
with energy or primary particle type was seen. When using QGSJetII 
or Sibyll2.1, Corsika and Conex predict the (X max ) with a difference 
smaller than 7 g/cm 2 , and the RMS(X max ) with a difference smaller than 
5 g/cm 2 . The differences in the slopes of a linear fit to the evolution of 
the (X max ) and RMS(X max ) with energy for Corsika and Conex are quite 
small (< 3 %). 

No assumption is made here for the cause of these differences. An in- 
vestigation for the possible cause could be done, but in the meanwhile these 
differences between the programs should be considered as systematics error 
in the analysis of the (X max ) and RMS(X max ) when one tries to infer the com- 
position abundance. Table [3] shows the suggested systematic uncertainties 
for the model predictions. Maximum values for the systematic uncertainties 



can be extracted from figures 4e, 5e, 6e and 7e 



Hadronic Model 


Mass (A) 


systematic uncertainty suggested for the model predictions 


(-^max) 

g/cm 2 


RMS(X max ) 
g/cm 2 


Elongation rate 
g/cm 2 per energy decade 


Sibyll2.1 


1 


2.57 - 1.05 X log 10 (E/EeV) 


-4.58 - 0.66 x log w {E/EeV) 


< 1.30 


55 


2.57 - 0.09 x log 10 (E/EeV) 


-7.78 - 0.52 x log 10 (E/EeV) 


< 0.08 


QGSJetII 


1 


3.73 - 0.31 x log w (E/EeV) 


-3.82 - 0.29 X log u) (E / EeV) 


< 0.20 


55 


5.13 - 0.70 x log 10 (E/EeV) 


-7.38 - 0.30 x log w (E/EeV) 


< 0.60 



Table 3: Systematic uncertainties suggested for the model predictions. 



Section [4] shows the parametrization of the X max distributions as a func- 
tion of energy and mass. The curves shown there can be used to estimate 
the first and second moments of the X max distribution from abundance cal- 
culations based on astrophysical arguments. As an example of the usage of 
this parametrization we have taken the astrophysical models developed by 
Berezinsky et al. [?] (Model 3) and Allard et al. [?] (Model A) and used our 
paremetrization to transform the abundance curves predicted by the models 
into a X max distribution. Figure 15 shows a X max distributions predicted by 
the models in comparison to the data measured by the Pierre Auger Observa- 
tory [?]. We have convolved the model predictions with a Gaussian detector 
resolution of 20 g/cm 2 . The utility of the parametrization is such that the 
models can be compared to the X max distribution instead of only the (X max ) 
and RMS(X max ). 
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(a) THIN = 1CT 5 (b) THIN = 1(T 6 (c) THIN = 1CT 7 

Figure 1: Longitudinal development of proton showers with energy 10 19 eV 
for different thinning factors. Each figure shows the development of one 
shower. Sib YLL2.1 was used for this study. 
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Figure 2: X max distribution for 1000 proton showers with 10 19 eV simulated 
with different thinning factors. Sibyll2.1 was used for this study. 
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(a) Conex- (X max ) versus energy (b) Corsika- (X max ) versus energy 
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(c) Conex- (X max ) versus mass, (d) CORSIKA- (X max ) versus mass. 
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Figure 4: (X max ) as a function of energy and mass as calculated by Corsika 
and Conex using Sibyll2.1. Showers have been simulated with primary 
energy ranging from 10 17 to 10 20 ' 4 eV in steps of log 10 (.E/e^) = 0.1 and 
primary nuclei types with mass: 1, 5, 15, 25, 35, 45 and 55. A set of 1000 
showers has been simulated for each combination. Not all energies and pri- 
maries are shown for clarity. 
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(c) Conex- (X max ) versus mass, (d) Corsika- (X max ) versus mass. 







o A 
A 


= 1 

— 15 

= 55 1 

f 10 
; 5 

I 

i 

x -5 
-10 




fllll 




rti 






1 1 


I'' 





17 17.5 18 18.5 19 19.5 20 20.5 15 0.2 0.4 0.6 0.8 1 1.2 1.4 tJB 1.6 

log (E/eV) I0 9, < A > 



(e) Difference between CORSiKA(f ) Difference between CORSIKA and 
and Conex versus energy. Lines areCONEX versus mass, 
linear fit to lp and 55p data. 

Figure 5: (X max ) as a function of energy and mass as calculated by Corsika 
and Conex using QGSJetII. Showers have been simulated with primary 
energy ranging from 10 170 to 10 20 ' 4 eV in steps of \og 1Q (E/eV) = 0.1 and 
primary nuclei types with mass: 1, 5, 15, 25, 35, 45 and 55. A set of 1000 
showers has been simulated for each combination. Not all energies and pri- 
maries are shown for clarity. 
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(a) Conex- RMS(A" max ) versus en-(b) Corsika- RMS(A" max ) versus 
ergy. energy. 
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(c) Conex- RMS(X n 
mass. 
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versus(d) Corsika- RMS(A" max ) versus 
mass. 




(e) Difference between Corsika(i ) Difference between CORSIKA and 
and Conex versus energy. Lines areCONEX versus mass, 
linear fit to lp and 55p data. 

Figure 6: RMS(X max ) as a function of energy and mass as calculated by 
Corsika and Conex using Sibyll2.1. Showers have been simulated with 
primary energy ranging from 10 17 ' to 10 20 ' 4 eV in steps of log 10 (E/eV) = 0.1 
and primary nuclei types with mass: 1, 5, 15, 25, 35, 45 and 55. A set of 
1000 showers has been simulated for each combination. Not all energies and 
primaries are shown for clarity. 
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(a) Conex- RMS(X max ) versus en-(b) Corsika- RMS(X max ) versus 
ergy. energy. 
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versus(d) Corsika- RMS(A" max ) versus 
mass. 



— I5r 

1 : 
S n~- 

2 : 
I 

8 » «7 



rf t j"'P' t '{| l M' 1 W^ u 



19 19.5 20 20.5 
log (E/eV) 




(e) Difference between Corsika(i ) Difference between CORSIKA and 
and Conex versus energy. Lines areCONEX versus mass, 
linear fit to lp and 55p data. 

Figure 7: RMS(X max ) as a function of energy and mass as calculated by 
Corsika and Conex using QGSJetII. Showers have been simulated with 
primary energy ranging from 10 17 ' to 10 20 ' 4 eV in steps of log 10 (E/eV) = 0.1 
and primary nuclei types with mass: 1, 5, 15, 25, 35, 45 and 55. A set of 
1000 showers has been simulated for each combination. Not all energies and 
primaries are shown for clarity. 
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Figure 8: RMS(X max ) versus (X max ) for all primary particles used in this 
work. The corresponding RMS(X max ) and (X max ) for a nuclei with mass 
55 has been taken as reference. Each blob corresponds to the (X max ) and 
RMS(X max ) predictions for one primary particle at different energies. 
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(a) Slope of (X max ) versus mass -(b) Slope of (X max ) versus mass 
SIBYLL2.1. QGSJetII. 
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(c) Slope of RMS(X max ) versus(d) Slope of RMS(X max ) versus 
mass - Sibyll2.1. mass - QGSJetII. 

Figure 9: Slope Analysis. Slope of a straight line fit to the variation of (X max ) 
and RMS(X max ) with energy (plotted as a function of mass). Slope in units 
of g/cm 2 /eV. 
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Figure 10: Parametrization of the X max distribution. CONEX- SIBYLL2.1. 
These figures show the general behavior of the three parameters used to 
describe the X max distributions as a function of energy and mass. 
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Figure 11: Differences in (X max ) versus energy. Comparison between the 
simulation, the direct fit of the X max distribution using equation [2] and the 
calculation using equation [3] and table [T| 
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(a) Conex- Sibyll2.1 



(b) Conex- QGSJetII 
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Figure 12: Differences in RMS(X max ) versus energy. Comparison between 
the simulation, the direct fit of the X max distribution using equation [2] and 
the calculation using equation [3] and table [T] 
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Figure 13: Comparison of (X max ) and RMS(X max ) for Corsika and Conex 
function of energy. 
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(a) A = 1 - 10 18 eV. 




(b) A = 1 - 10 19 eV 



(c) A = 1 - 10 2U eV. 
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(d) A = 55 - 10 18 eV. (e) A = 55 - 10 19 eV. (f) A = 55 - 10 2U eV. 

Figure 14: Simulated X max distribution fitted by a Gaussian convoluted with 
an exponential (equation [2j. In this example we show showers simulated with 
Conex and Sibyll2.1. Full line shows the X mgx distribution and dashed 
line the fit of equation |2l 
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(c) Allard et al. - CORSiKA-(d) Berezinsky et al. - CORSIKA- 
QGSJetII. QGSJetII. 

Figure 15: X max distributions. Data measured by the Pierre Auger Observa- 
tory with energy 10 18 ' < E < 10 18,1 eV [?]. Astrophysical models extracted 
from [?, ?]. The models have been calculated at E — 10 18 - 05 eV. The curves 
have been calculated using the parametrizations proposed above. 
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